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ABSTRACT 

We present the results of Smoothed Particle Hydrodynamics (SPH) simulations of 
the formation of a massive counterrotating disk in a spiral galaxy. The current study 
revisits and extends (with SPH) previous work carried out with sticky particle gas dy- 
namics, in which adiabatic gas infall and a retrograde gas-rich dwarf merger were tested 
as the two most likely processes for producing such a counterrotating disk. We report 
on experiments with a cold primary similar to our Galaxy, as well as a hot, compact 
primary modeled after NGC 4138. We have also conducted numerical experiments with 
varying amounts of prograde gas in the primary disk, and an alternative infall model 
(a spherical shell with retrograde angular momentum). The structure of the resulting 
counterrotating disks is dramatically different with SPH. The disks we produce are con- 
siderably thinner than the primary disks and those produced with sticky particles. The 
time-scales for counterrotating disk formation are shorter with SPH because the gas 
loses kinetic energy and angular momentum more rapidly. Spiral structure is evident in 
most of the disks, but an exponential radial profile is not a natural byproduct of these 
processes. The infalling gas shells that we tested produce counterrotating bulges and 
rings rather than disks. The presence of a considerable amount of preexisting prograde 
gas in the primary causes, at least in the absence of star formation, a rapid inflow 
of gas to the center and a subsequent hole in the counterrotating disk. For a normal 
counterrotating disk to form, there must either be little or no preexisting prograde gas 
in the primary, or its dissipative influence must be offset by significant star formation 
activity. The latter scenario, along with the associated feedback to the ISM, may be 
necessary to produce a counterrotating disk similar in scale length and scale height to 
the primary disk. In general, our SPH experiments yield stronger evidence to suggest 
that the accretion of massive counterrotating disks drives the evolution of the host 
galaxies towards earlier (SO/Sa) Hubble types. 

Subject headings: galaxies: spiral — galaxies: structure — galaxies: evolution — galax- 
ies: interactions — galaxies: kinematics and dynamics — hydrodynamics 
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1. Introduction 

There are only a handful of known cases of 
massive counterrotating disks in spiral galaxies 
to date, and yet counterrotation in spirals can- 
not be deemed a rare phenomenon by any stan- 
dards. There are hints that it may be quite 
common, in fact, in early-type spirals, partic- 



ularly SOs (Kuijken, Fisher & Merrifield 1996). 
The origin of any counterrotating mass (gas or 
stars) within a spiral disk is an important un- 
solved problem with profound implications for 
the formation and evolution of all spiral galaxies, 
but the existence of a significant retrograde mass 
component (comprising anywhere from ~ 10-50% 
of the total mass of the disk system) is a partic- 
ularly intriguing question that threatens to rad- 
ically alter our view of the evolution of spiral 
galaxies. The rogues' gallery of spirals with such 
massive counterrotating disks currently boasts 



the following members: NGC 4550 ([Rubin, Gra- 
[ham fc Kenney 1992| ), NGC 7217 (jvierrifield &[ 



[uijken 1994[), N GC 4826 ( [Braun, Walterbos k 



[ennicutt 199^ ), NGC 3626 ([Ciri, Bettoni fc 



Palletta 1995[) , NGC 3593 ([Bertola et al. 1996[) , 
and NGC 4138 ( [Tore, Broeils fc Haynes 1996[) . 

Apart from the challenge they present to the 
traditional view of spiral galaxies that has evolved 
over the last few decades, counterrotating disks 
raise several questions about other astrophysi- 
cal processes, such as the role and fate of gas 
in galaxy interactions, star formation in galax- 
ies that contain counterrotating gas, the accre- 
tion rates and star formation histories of spi- 
ral galaxies in general, and the impact of coun- 
terrotating populations on the overall stability 
of the disk system. Even though a recent sur- 



vey of SOs ( Kuijken, Fisher &: Merrifield 1996 ) 
found counterrotating gas in almost a quarter of 
the sample, none of these galaxies have counter- 
rotj ating stars. — Why is tho counterrotating gas 



not ! forming stars? If star formation is inhibited 
in counterrotating disks, how does one explain 
NGC 4550 and others with stellar counterrotat- 
ing disks? 



uct of the galaxy formation process. The theory 
of formation of a spiral galaxy from a spinning 
protogalactic cloud does not admit the possibility 
of bidirectional spin being imparted to the disk 
system. Subsequent accretion or merger events 
are a much more plausible explanation, and even 
these are severely constrained by the observed 
coldness of the counterrotating galaxies. Dissi- 
pationless mergers, especially between progeni- 
tors with comparable masses, can be ruled out 



in the general case, although [Pfenniger (1998) 



has recently been able to produce a remnant re- 
sembling NGC 4550 with a collisionless merger of 
two spirals with special initial conditions. This 
leaves minor, gas-rich mergers, and gas accretion 
or infall, as the most promising candidates. 

The most puzzling aspect of massive counter- 
rotating disks in spiral galaxies is that the host 
galaxies appear quite normal in every other re- 
spect, and there is no evidence of excessive thick- 
ening of the primordial disks due to the accre- 
tion of the counterrotating disk. This suggests 
that the accretion process must not be a rapid 
or violent one, but there may be deeper implica- 
tions here for the interaction histories of all spiral 
galaxies. A recent study of the tidal thickening of 
galaxy disks (Reshetnikov & Combes 1996) indi- 
cates that the ratio of scale length to scale height, 
h/zo, is 1.5-2 times lower for interacting disks. 
But if any galaxy can be assumed to have under- 
gone an interaction in the past, the thinness of 
the non-interacting sample proves that this ratio 
returns to its higher value after a certain amount 
of time (order of 1 Gyr). If it can be proved 
that a spiral can double its disk mass without de- 
stroying the primordial disk in the process, then 
the present-day appearance of spiral disks can no 
longer preclude such interactions in their past. 

The claim that the thinness of spiral disks 
places a stringent limit on past accretion ( p?6th 
fc Ostriker 199^ ) is further challenged by re- 



L is very unlikely LhaL couiiLerroLaLiiig sys- 



cent simulations showing that the halo absorbs 
a good portion of the orbital energy and angular 
momentum of the satellite in spiral-dwarf merg- 
ers ([Walker, Mihos Sz Hernquist 1996 ; Huang & 



terns can be produced mdigenously or as a byprod- 



Carlberg 1997), and analytical results indicating 



that an isothermal halo may even shield the disk 
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from an external tidal field due to a satellite ( Mu- 
ral _&_TremBine_J_997). A study of the lopsid- 
edness of the disks of field spirals suggests that 
the accretion rate of spiral galaxies may be as 
high as one small (~10% mass) companion ev- 
ery 4 Gyr ( [Zaritsky fc Rix 1997| ). This is stih 
a very uncertain estimate, and a more accurate 
estimate of the accretion rate, and a knowledge 
of how much accretion a disk galaxy can with- 
stand, are important questions from a cosmolog- 



ica t point of view. Observations of counterrotat- 



ing disks combined with a good understanding 
of how they form provide a striking new way to 
independently constrain these estimates, because 
the accreted matter can be easily distinguished 
in a counterrotating galaxy. 

With the aim of understanding the origin of 
massive counterrotating disks in spirals and SOs, 
we have developed a numerical code and run 
hydrodynamical simulations to investigate the 
processes that are most likely to produce such 
bizarre systems. We have combined an N-body 
gravity solver with a gas dynamics particle code 
for this purpose. To study the basic parame- 
ters of the processes involved, we first adopted 
a quick-and-dirty "sticky particle" gas dynam- 
ics approach that allowed us to test various sce- 
narios with relatively small investments of CPU 
time. In this way we were able to test massive 
counterrotating disk formation in a fiducial cold 
primary (IThakar fc Ryden 1996| , hereafter TR) 
as well as model the formation of the recently 
discovered counterrotating disk in the early-type 
spiral NGC 4138 ([Thakar et al. 1997| , hereafter 
TRJB). For each type of primary modeled, we 
have tested two theories of origin: adiabatic (sec- 
ular) gas infall and a gas-rich dwarf merger. In 
order to produce a disk with opposite spin, both 
of these processes require a retrograde orbit of 
accretion for the infalling gas or satellite galaxy 
with respect to the primary's spin. Gas infall 
works well for both late and early-type primaries, 
but a dwarf merger, especially with a substantial 
amount of dissipationless matter in the dwarf, is 
not viable for a cold primary because it plays 
havoc with the primary's disk. 

An accurate representation of the gas in astro- 



physical systems typically requires that the hy- 
drodynamical conservation equations be solved 
and the effects of physical processes such as 
shocks and viscosity be included. This requires 
the gas to be modeled as a fluid, but a good com- 
promise can be achieved with a particulate rep- 
resentation if each gas particle is smeared over a 
finite volume and its physical properties averaged 
or smoothed over that volume. This Smoothed 
Particle Hydrodynamics (SPH) appr oach ([Lucy 
19771 ; iGingold k, Monaghan 1977] ; [Monaghan 



1992) meshes well with an N-body particle code. 



but is much more costly in terms of CPU re- 
sources than the sticky particle approach. We 
therefore reserved it for a more detailed and re- 
stricted look at the structure of the counterro- 
tating disk formed by reexamining (with SPH) 
a carefully chosen subset of the simulations we 
presented in TR and TRJB. We present the re- 
sults of these simulations here. Although most 
of the simulations are reruns of those presented 
before with sticky particles, a couple of new runs 
are also discussed. Our SPH code currently does 
not include star formation. We intend to incor- 
porate this in our simulations in the future, but 
as we discuss in §4 below, we do not expect it to 
have a profound impact on the structure of the 
counterrotating disks formed. 

Our computational method is described in the 
following section. Results of simulations are pre- 
sented in §3, followed by discussion and conclu- 
sions in §4 and §5 respectively. 

2. Computational Method 

Our SPH code is an implementation in the C 



language of the TREESPH code (Hernquist & 



Katz (1989)1 , hereafter HK). We started with a 



basic C version of the tree code (Barnes & Hut 



1986 ) that performs the force calculation, and 



optimized it for execution on the Cray Y-MP at 
the Ohio Supercomputer Center (see TR for de- 
tails). These changes included rewriting the tree 
construction and traversal routines in such a way 
that they made extensive use of the vector ca- 



pabilities of the Y-MP (iHernquist 1990| ; [Makino 



1990 ). The outcome was to speed up the execu- 
tion by a factor of 5 or more. We then added 
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the SPH code to it according to the prescription 



Lattanzio 1985, eq. 2.8 in HK). However, we 



given by HK, using the same vectorization tech- 
niques as before to optimize the tree traversals 
required for nearest neighbor searches. Vector- 
ization of the SPH calculations was tricky, and 
we benefited considerably from a perusal of the 
vectorized Fortran TREESPH code (a copy was 
kindly provided by Lars Hernquist). To test the 
final result, we ran the one-dimensional and 3- 
d tests described in HK to ensure that it gave 
identical results. 

For the simulations described here, we have 
assumed an isothermal gas with a cutoff temper- 
ature of 10^ K. This has been shown to be a fairly 
good assumption for the ISM in disk galaxies, 
where most of the gas hovers close to the cutoff 
temperature due to the short time-scale of the 
radia tive cooling process ([Barnes &: Hernquist 



19^; iMihos fc Hernquist 1994a| , 1994c). The mo- 
tivation for assuming an isothermal gas is that it 
yields a tremendous saving in the CPU resources 



do not use it for gravity softening. Our grav- 
ity solver uses the standard softening scheme in 
which the gravitational potential is of the form 
$ oc (r^ -|- e^)~^/^ (corresponding to a Plummer 
density profile), with softening length e = 0.5 kpc 
for all simulations. Although our code provides 
the capability to use different softening lengths 
for the collisionless and gas particles, we have 
used the same length for both species in the sim- 
ulations reported here. The softening length ef- 
fectively limits our resolution, but it is also nec- 
essary to suppress the two-body effects in our 
moderate- simulations. 

Our time-integration scheme is different from 
the one used by HK. Their scheme requires man- 
ual synchronization of the particle positions and 
velocities by applying midpoint corrections (eqs. 
2.42 - 2.45 of HK). We decided instead to adopt 
a slightly different scheme that does not re- 
quire these corrections ( [Katz, Weinberg fc Hern" 



rec juired for each simulation, since solving the en- quist 199^ ). Although this scheme is not time- 



ergy equation is a non- trivial exercise that re- 
quires a semi-implicit scheme (see HK) and must 
be performed at the smallest time-scale involved 
in the simulation. 

Even with the isothermal assumption, and 
with the maximum possible vectorization, our 
SPH code typically requires 5-10 times the num- 
ber of CPU resources required for the sticky par- 
ticle code. This is mostly because the gas par- 
ticles have their own individual time-steps (each 
particle's time-step is a power-of-two subdivision 
of the system time-step, so that all gas particles 
are synchronized at the end of a system time-step 
as discussed in HK), and this causes the perfor- 
mance to degrade when there are even a few par- 
ticles in the lowest time-step bin. To keep this 
problem under control, we do impose a minimum 
smoothing length (usually about a quarter of the 
softening length), but even so the isothermal as- 
sumption is more of a necessity than an option 
for us. As we discuss below, this assumption 
does present some problems when we attempt to 
model a gas-rich dwarf galaxy. 

We have adopted the spherically symmetric 



spline kernel for SPH averaging ( Monaghan fc 



symmetric result of the particle positions 
and velocities being advanced in a different way 
than in HK, it is simpler and faster without in- 
troducing any additional error. 

For all simulations reported here, we use the 
combined form of the bulk and Von Neuman- 
Richtmyer artificial viscosity given by eqs. 2.24- 
2.25 in HK, and applied the restriction necessary 
to avoid viscous cooling of receding particles as 
specified by HK. This viscosity introduces less 
shear than other artificial viscosities, although it 
does not provide as accurate a description of the 
flow near shocks. We felt that it was more im- 
portant to avoid excessive shear in disk galaxy 
simulations, and therefore adopted the combined 
form of viscosity. 

We use 32k particles each to represent the pri- 
mary's halo and disk, and up to 20k particles for 
the gas depending on its mass. Where a dwarf 
galaxy is involved, the stars and dark matter in 
the dwarf are represented by enough particles so 
that the per particle mass is roughly the same as 
in the case of the primary. A nagging problem 
with disk simulations is the disk heating due to 
the graininess of the potential that is inherent 
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in a particle code (Hockney & Eastwood 1989). 



The problem is more severe in our simulations be- 
cause of the longer time-scales involved in retro- 
grade interactions. Even with a softening length 
as high as 0.5 kpc, we find that there is consider- 
able disk heating over the duration of our simu- 
lations with the number of particles we have cho- 
sen, and the heating of the disk due to the inter- 
action itself is always of the order of the numer- 
ical heating. Unfortunately, the two are not sep- 
arable, and this limits our ability to estimate the 
disk thickening accurately. Although we would 
like to use more than 32k particles each to rep- 
resent the disk and halo, our CPU resource lim- 
itations currently prevent us from running con- 
siderably larger N simulations. 

The primary galaxy physical parameters are 
the same as we used before for the cold primary 
(TR) as well as for NGC 4138, the hot primary 
(TRJB). For convenience, they are repeated in 
Table |. Our fiducial cold primary parameters 
are typical for an Sbc spiral like the Milky Way, 
and the disk is bar-unstable as is consistent with 
its coldness (Binney &: Tremaine 1987). Unless 



explicitly stated, the primary disk in most of the 
simulations reported here contains no primordial 
gas. The halo is a truncated isothermal sphere 
(truncated at r = R^) with a logarithmic profile 
and core radius Rc as per eq. 2-54a in Binney 
& Tremaine (1987). Although most of the SPH 
simulations are reruns of previous sticky particle 
runs, there are some changes and some new sim- 
ulations. The input parameters for the gas infall 
and merger simulations are listed in Tables § and 
^ respectively. Correspondences to the original 
simulations in TR and TRJB are pointed out, 
where applicable, when we discuss the results in 
the following section. To avoid confusion with 
the sticky particle runs, we have named the SPH 
gas infall runs G1-G8, and the SPH dwarf merger 
runs Dl and D2. 

The models for gas infall and mergers are the 
same as in TRJB. The infalling gas is configured 
as a long rectangular column of uniform den- 
sity and square cross-section. The gas columns 
we use in the SPH simulations may be shorter 
than the corresponding sticky particle simula- 



tions, primarily because the SPH runs require 
much longer to run and hence a shorter col- 
umn with a correspondingly shorter infall time 
is preferable. The primary disk plane is selected 
to be the xy-plane, and the long axis of the in- 
falling gas slab is parallel to the y-axis. The pri- 
mary disk rotation is counter-clockwise, and the 
center of the primary disk is at the origin of the 
coordinate system. The gas is then placed in the 
first quadrant in the same plane as the primary 
disk, and is given an initial kick mostly in the 
negative y direction (an initial speed in the posi- 
tive X direction is also given in order to increase 
the initial angular momentum of the gas in some 
cases, see Table Q). This ensures that the gas is 
accreted on a retrograde orbit with respect to the 
primary disk spin. 

For the dwarf mergers, the setup is basically 
the same except that the dwarf galaxy replaces 
the gas column in the first quadrant. The dwarf 
galaxy model is the same as described in TR 
and TRJB. We have only attempted two gas-rich 
dwarf mergers (Dl and D2), both with NGC 4138 
as the primary, since they were not feasible for 
the cold primary due to the problems discussed 
in TR. Both mergers have identical input param- 
eters and are a repetition of merger M2 from 
TRJB. SPH is used from the beginning for Dl, 
whereas for D2 the merger is run with sticky par- 
ticles until the dwarf galaxy is completely dis- 
rupted by the primary, at which point SPH re- 
places sticky particles for further calculations. 
This is due to the difficulties we encountered in 
using isothermal SPH to model the dwarf galaxy 
(see §3.2.2). 

3. Results 

The primary galaxy is evolved in isolation for 
a few dynamical times before the gas or dwarf 
galaxy is introduced. Hence the simulation re- 
sults are shown starting ai t 1.5 Gyr for the 
cold primary, and t ;^ 0.5 Gyr for NGC 4138 
(which has a smaller dynamical time). The sim- 
ulation stops once most of the gas has settled in 
the plane of the primary and the evolution of the 
gas disk has tapered off to nearly zero. 
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Table 1: Parameters for Primary Galaxy. 







Disk 






Halo 






Mass 


Radius 


Scale Len. 


Scale Ht. 


Mass 


Radius 


Core Rad. 




Md [Mq] 


Rd [kpc] 


Rd [kpc] 


hd [kpc] 


Mh [Mq] 


Rh [kpc] 


Rc [kpc] 


Cold Primary [TR] 


5.5 X 10^" 


21.0 


3.50 


0.3 


2.2 X 10^^ 


50.0 


1.0 


NGC 4138 [TRJB] 


2.0 X 10^° 


6.0 


1.25 


0.2 


1.5 X 10^1 


20.0 


5.0 



Table 2: Parameters for Cas Infall Simulations. 





M/ 


T/ b 






^0^ 






Gl 


4.40 


0.3 


0.10 


50 


50 


300 


20 


G2 


4.40 


0.3 


0.00 


50 


50 


300 


10 


G3 


1.32 


0.3 


0.05 


25 


25 


150 


10 


G4 


1.32 


0.3 


0.05 


25 


5 


100 


10 


G5 


1.45^ 


0.3 


0.10 


25 


10 


100 


10 


G6 


1.76t 


0.3 


0.05 


25 


5 


100 


10 


G7 


1.06^ 


0.2 






30 




10 


G8 


1.06^ 


0.1 






30 




10 



''total gas mass [10^° Mq] 

''initial velocity of the gas in the negative y direc- 
ti on, specified in units of the centripetal velocity 
yjG{Mh + Md)/Xo, where Mh + Md is the total mass 
of the primary (halo+disk) 

'^initial velocity of the gas in the positive x direction, spec- 
ified as a fraction of Voy 

''initial y-distance of the gas column from the center of the 
primary disk [kpc] 

"initial a;-distance of the gas column from the center of the 
primary disk [kpc] 

'length of the infalling gas column [kpc] 
^thickness of the gas column (side of square cross section) 
[kpc] 

^G5 includes a prograde ring of gas in the primary with 
~ 10% of the mass of the infalling gas (hence Mg = 
1.32 + 0.13), and G6 includes a prograde disk of gas in 
the primary with ~ 20% of the mass of the stellar pri- 
mary disk (hence Mg = 1.32 -|- 0.44) 

*G7 and G8 are infalling gas shells having retrograde angu- 
lar momentum with rotational speed Voy, radius Xq and 
thickness Tq. 

3.1. Fiducial Cold Primary 

For the cold primary, our first infall simula- 
tion is a repetition of the second continuous in- 
fall simulation discussed in TR, but in the SPH 



Table 3: Dwarf Galaxy Parameters for Mergers. 

Af/ M,^ Mp' F/ R,' R/ Rh^ 

[10'"Mq] [kpc] 

D1/D2 1.06 0.26 2.64 0.27 4.0 4.0 8.0 



''mass of gas 

''mass of stars 

'^mass of dark matter 

''ratio of gas mass to total mass 

''radius of stellar sphere 

'radius of gas sphere 

^radius of halo 

run we use a shorter column of gas. The infall 
simulation shown in Fig. |l| has all other inputs 
the same as before (see Table |2| for input param- 
eter values). A counterrotating disk with spiral 
structure has started to form by i ~ 4.5 Gyr. The 
side view confirms the disk formation and indi- 
cates that the counterrotating disk is uniformly 
thin throughout. Since the width of the last two 
panels in Fig. |l| is 40 kpc, the size of the counter- 
rotating disk appears to be comparable to that 
of the primary (radius 21 kpc). One has to be 
careful in making this comparison, however, since 
the radial profile of the counterrotating disk is 
not close to that of the primary. Saying that the 
sizes are comparable does not mean that the scale 
lengths are comparable. The lack of a centrally 
concentrated or exponential radial profile is clear 
from both the side and top views for t '^b Gyr. 
There is a hole in the inner regions of the coun- 
terrotating disk, and this, along with the larger 
size, could be the result of the gas having too 
much angular momentum initially. 

Although the observed instances of counter- 
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rotating disks do not show a consistent radial 
profile, it is nonetheless important to determine 
whether a counterrotating disk with a similar 
scale-length and profile as the primary disk can 
be formed. There is at least one example of this, 
NGC 4550 ( [Rix et al. 1992^ . We have exper- 
imented with lower values of the initial angu- 
lar momentum of the infalling gas to determine 
what effect it has on the radial profile and struc- 
ture of the counterrotating disk formed. Fig. ^ 
presents the results of run G2, which has lower 
initial angular momentum of the gas as well as a 
narrower (and hence denser) column of infalling 
gas than Gl. Since the panel widths in Fig. ^ are 
the same as in Fig. |^, the differences are easy to 
see. The size and structure of the counterrotat- 
ing disk formed are quite different. The disk is 
smaller (by about 25-30%) than the disk formed 
in Gl, and the radial profile is more centrally 
concentrated. There is a lot more structure visi- 
ble as the gas has piled up in several places due 
to shocks. The side view shows that the disk is 
as thin as in Gl, and it is significantly thinner 
than the primary disk. 

A comparison of the primary disk at the be- 
ginning and the end of G2 with the counterro- 
tating disk is shown in Fig. ^. The primary disk 
shows a bar even before the gas is introduced, 
but the bar does not become more pronounced 
after the formation of the counterrotating disk. 
It is, in fact, obscured by the inclination of the 
disk in the top view of the final timestep. The 
final thickness of the primary disk is hard to mea- 
sure because of its inclination, but we do measure 
it by compensating for the inclination and com- 
puting the z = plane at each location in the 
disk (see below). However, it is obvious from 
Fig. I that the primary disk thickness has more 
than doubled, and it now resembles a lenticular 
(SO) disk rather than an Sb disk. This is further 
corroborated by the rotation curves and velocity 
dispersion plots for the final disks shown in §4, 
where we also discuss some of the possible rea- 
sons for the pronounced heating experienced by 
the cold primary. 

The counterrotating disk is of comparable size 
but is less than half as thick as the primary disk. 



The velocity fields of the primary disk and the 
counterrotating disk are shown in Fig. ^. The pri- 
mary's velocity field is considerably hotter than 
the velocity field of the counterrotating disk, but 
they are clearly antiparallel. 

The radial profiles of the primary disk and 
the counterrotating disks at the end of simula- 
tions Gl and G2 are compared with each other 
in Fig. 1^. The dip in the particle density between 
~ 1 — 3 kpc is clearly visible for Gl, and although 
there is a slight dip at those radii even for G2, the 
G2 profile comes much closer to the exponential 
radial profile of the primary disk. 

The thickness of the primary disk is plotted for 
simulations Gl and G2 in Fig. |6[ A comparison 
of the thickness plots with those in TR shows 
some similarities. The thickness of the post- 
accretion primary disk is only slightly higher 
than the thickness of the isolated primary disk 
evolved over the same time period. The actual 
thicknesses for Gl and G2, however, are signifi- 
cantly higher (by about 25%) than the thickness 
plotted for continuous infall in TR. This is prob- 
ably due to the fact that we are using a smaller 
softening length (0.5 kpc instead of 1.0 kpc used 
by TR), even though we have doubled the num- 
ber of particles used to represent the halo for 
the simulations in this paper (we used 16k par- 
ticles for the halo in TR). The thickness curve 
for the isolated disk is also higher by a similar 
amount compared to the corresponding curve in 
TR. Some of the increase in the heating of the 
primary can also be attributed to the fact that 
the infalling gas column is shorter and hence 
more dense than the one used in TR, and the 
time-scale for formation of the counterrotating 
disk has decreased by ^ 30% (see §4.1). This 
gives the primary disk less time to adjust. Fur- 
thermore, the smaller softening length aggravates 
the gravitational impact of the incoming gas on 
the primary. 

We are not so concerned with the heating of 
the primary in the simulations reported here, be- 
cause we have established in TR and TRJB that 
if the rate of accretion is low enough, the effect on 
the primary can be kept to an acceptable level. 
As mentioned above, our motivation for hasten- 
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ing the counterrotating disk formation process 
by using a shorter column of gas than before was 
to save CPU time expended on each simulation. 
The fact that this causes considerably more heat- 
ing of the primary underlines the sensitivity of 
the cold primary to the rate of accretion. 

3.2. Hot Primary - NGC 4138 

3.2.1. Gas Infall 

The first infall simulation for NGC 4138, 03, 
is a repetition of simulation G4 from TRJB, al- 
though the total mass of the gas is 25% higher 
in the SPH runs. The remaining parameters are 
the same as in the sticky particle run. The re- 
sults are shown in Fig. |^. The counterrotating 
disk has started to form at t ~ 1.8 Gyr, and af- 
ter t ~ 2.6 Gyr, it does not evolve very much. 
The size of the counterrotating disk is signifi- 
cantly larger than the primary since it extends 
over more than two-thirds of the width of the 
panel (20 kpc) in the side-view at t > 2.6 Gyr, 
but it clearly does not have an exponential radial 
profile. There is a small central mass concentra- 
tion. The lack of an exponential radial profile is 
easier to see in the side view at t = 3.0 Gyr, be- 
cause the disk is slightly tilted to the line of sight, 
and an outer ring is visible in addition to the 
central concentration, with a drop in the particle 
density at intermediate radii. The radial density 
distribution is plotted below in Fig. 10. As in 



the case of the cold primary (Gl), the lack of an 
exponential profile indicates that a good fraction 
of the particles do not lose enough angular mo- 
mentum to allow them to settle in the inner few 
kpcs of the disk. Spiral structure is evident but 
it is not very strong. The side view at t = 2.6 
Gyr shows the counterrotating disk to be quite 
thin, although this is hard to see in the last panel 
because of the inclination of the disk. 

The velocity fields of the primary and counter- 
rotating disk are compared in Fig. ^. The larger 
size of the counterrotating disk is also evident in 
the velocity field plot. The velocity field of the 
counterrotating disk is colder than that of the 
primary and more regular. 

A lower angular momentum version of G3 is 



attempted in G4, shown in Fig. ^. The initial 
velocity of the gas column is the same as G3, 
but the initial distance of the gas from the cen- 
ter of the primary is smaller by a factor of 5 
(see Table §). The length of the gas column 
is also two-thirds of that in G3. The parallels 
between G1/G2 and G3/G4 are quite obvious. 
Once again, the lower angular momentum pro- 
duces a counterrotating disk with a more cen- 
trally concentrated radial profile and a smaller 
overall size. A comparison of the radial profiles 
of G3 and G4 with each other and the primary's 
radial profile is shown in Fig. 10. Even though 



the radial profile of the counterrotating disk in 
G4 is not very close to the exponential profile of 
the primary, it is much closer to that than G3. 
There is no pronounced dip in the particle den- 
sity in the inner radii for G4, as there is for G3. 

There is not much difference in the heat- 
ing (and thickening) of the primary disk with 
the SPH simulations of NGC 4138 compared to 
TRJB. Although the heating is a little higher on 
average in our SPH simulations, this is due to the 
fact that we are using shorter gas columns. This 
makes the accretion rate higher. In the cases 
where a lot of gas does make it to the nuclear 
portions of the primary, there is slightly more 
heating of the primary. In the other cases, the 
levels are comparable to those obtained in TRJB 
(after allowing for the shorter gas columns), and 
hence the thickness plots are not repeated here. 

3.2.2. Gas-Rich Dwarf Merger 

Our previous studies have established that the 
gas-rich dwarf merger model of counterrotating 
disk formation is only feasible for a hot primary 
like NGC 4138. TR found that even a merger 
with a 10% mass satellite produced unacceptable 
levels of heating in a cold primary disk, whereas 
the hot primary used by TRJB fared much better 
in this regard and was able to withstand a gas- 
rich dwarf merger under restricted conditions. 

Apart from the issue of whether a dwarf merger 
can produce a counterrotating disk, we encoun- 
tered a technical difficulty with our gas-rich dwarf 
galaxy model. A stable model of a collisionally 
supported gas-rich dwarf galaxy is not easy to 
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obtain with an isothermal equation of state, be- 
cause the gas in the dwarf tends to shock and 
cohapse rapidly if the temperature is too low, 
and expand if it is too high. This inability of the 
isothermal model to yield a stable configuration 
is not very surprising, since isothermal pressure 
gradients do not grow rapidly enough to com- 



bat the gravitational collapse of the gas (Tohline 
19@. A rotationally supported dwarf would 
probably not alleviate this situation very much, 
because it would still collapse in the direction not 
affected by the centrifugal force, thereby yielding 
a quasi-equilibrium isothermal disk. As such, we 
restrict ourselves to testing only one model for 
which we have a reasonable expectation of suc- 
cess, and even for this we have to adjust the gas 
temperature to keep the dwarf galaxy in equilib- 
rium. 

We have selected merger M2 from TRJB, with 
a dwarf that is ~ 18% as massive as the pri- 
mary. Two versions of this merger are presented 
here. In the first case (run Dl), we use SPH 
with a higher isothermal temperature until the 
dwarf galaxy is tidally disrupted by the primary, 
and thereafter reduce the temperature of the gas 
to the usual value (10^ K). In the second ver- 
sion (run D2), we use the sticky particle results 
to evolve the dwarf until it forms a thick disk 
around the primary, and then switch to SPH with 
the usual settings. 

The results for Dl, with SPH being used all 
the way, are shown in Fig. 11. The top view 



shows that the dwarf makes one entire pass 
around the primary before getting completely 
stripped, and ends up forming a very barred 
counterrotating disk. In the side view the disk 
appears to be quite thin, and the velocity field 
diagrams shown in Fig. |l2| confirm that a disk is 
indeed formed. The stellar particles (not shown) 
of the dwarf form a thick flattened cloud around 
the primary, in the same way as they did with 
sticky particles. The velocity fields also indicate 
that the primary disk is heated up more than in 
the infall case (compare with Fig. 

When SPH is introduced at a late stage af- 
ter the dwarf has already been tidally stripped 
and has formed a thick disk around the pri- 



mary, the evolution of the disk is very different 
(Fig. |l3|). D2 dramatically illustrates the differ- 
ences between sticky particle gas dynamics and 
SPH. Within a short time 0.4 Gyr), the thick 
sticky particle disk has collapsed to a thin disk. 
There is a slow accumulation of mass in the cen- 
ter thereafter, but beyond the formation and dis- 
sipation of a few rings, the disk is fairly stable. 

It appears that the nature of the dissipation 
is very different with SPH. With sticky parti- 
cles, there is a slow but steady loss of kinetic en- 
ergy, whereas with SPH, shocking causes rapid 
loss of energy. The isothermal model proba- 
bly contributes to the speed of the energy loss, 
since there is no counteracting increase in ther- 
mal pressure to oppose the kinetic energy dissi- 
pation until the cutoff temperature is reached. 
Once the SPH gas forms a kinematically cold 
disk, evolution ceases as the dissipation drops 
off suddenly. The rapid loss of kinetic energy 
is demonstrated even more effectively when pro- 
grade gas is present in the primary (see below) 
and the counterrotating gas collides head-on with 
it. 

3.3. Other Experiments with NGC 4138 

3.3.1. Prograde Gas in the Primary 

Sticky particle simulation 16 (TRJB), which 
included a ring of prograde gas in the primary, 
showed that the prograde gas did not have much 
of an impact on the kinematics of the counterro- 
tating gas. Collisions between counterstreaming 
gas particles did send some gas to the center, 
but this was a modest effect and the prograde 
gas ring remained mostly intact. 



The SPH simulation G5 shown in Figs. 14 and 
15 is a repetition of 16 from TRJB. The average 
radius of the ring is 1.7 kpc and the width of the 
ring is 1 kpc (inner radius of ring = 1.2 kpc and 
outer radius = 2.2 kpc). The mass of the ring 
is ~ 10% of the mass of the infalling gas, and 
the ring has uniform density. The particles in 
the ring are given initial velocities appropriate for 
primary disk particles at those radii, and the ring 
is checked for stability for a couple of dynamical 
times prior to introducing the infalling gas. 
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We find that the behavior of the gas in the 
SPH version of this simulation is dramatically 
different from the sticky particle results. Vir- 
tually all of the prograde gas interacts with the 
counterrotating gas and falls to the center, leav- 
ing a gap in the counterrotating gas distribution 
at the original location of the prograde ring of 
gas. The interaction between the prograde and 
retrograde gas particles can be seen more clearly 
in the velocity fields plot. The angular momen- 
tum of the retrograde gas is reversed near the 
center of the primary as both the prograde and 
retrograde gas particles fall to the center. 

Going one step further, we included a pro- 
grade disk of gas in the primary along with the 
stellar disk. The mass of the gas disk is ~ 20% 
of the mass of the primary disk, and the radius 
of the gas disk is 4 kpc. As was the case for the 
prograde ring, the prograde disk is of uniform 
density, and we evolved it with the rest of the 
primary first for 200 Myr (without any counterro- 
tating gas) to ensure that it was stable on its own. 
Then we introduced the infalling gas at t = 0.7 



shown in Table ||. G8 has half the angular mo- 
mentum of G7. The results for G7 and G8 are 



Gyr. This is simulation G6 shown in Fig. 16. The 
behavior of the two gas components is similar to 
that seen in G5, only stronger. The prograde gas 
sweeps up all the retrograde gas that it comes 
in contact with, and both combine to form a nu- 
clear mass concentration within a Gyr or so. As 
in the case of the prograde ring, the velocity field 
plots (Fig. 17) show that the angular momen- 
tum of the infalling gas is reversed in the process. 



Ths rest of the infalling gas that does not come 



into contact with any of the prograde gas forms a 
counterrotating outer ring. The gap between the 
prograde nuclear disk/bulge and the outer coun- 
terrotating disk is larger in this case since the 
prograde gas occupied a larger area in G6 than 
in G5. 

3.3.2. Infalling Spherical Gas Shell 

Gas infall is not well-constrained by observa- 
tions, and as an alternative gas infall model, we 
tried a spherical shell of gas that is given ret- 
rograde angular momentum with respect to the 
primary disk. The two runs with such an infalling 
gas shell, G7 and G8, have the input parameters 



shown in Figs. 18 and 19. The shell collapses to 
a flat disk within a very short time (^ 0.5 Gyr). 
The speed of the collapse is somewhat higher for 
G8. The final size of the disk in G8 is about 
half the size of the disk in G7 (note that the 
last two panels of Fig. Il^ are half as wide as 



the corresponding panels in Fig. 18). The mass 
distribution of the two counterrotating disks is 
also different. There is an outer ring present in 
both disks, but only G8 shows a significant cen- 
tral mass concentration. The tilting of the disk 
at the end of G8 is in phase with the tilting of the 
primary disk. The formation of the outer ring in 
both cases, together with the dependence of the 
size of the ring on the initial angular momentum, 
suggests that the infalling shell mechanism is not 
particularly suitable for counterrotating disk for- 
mation. There is no a priori reason to believe 
that an infalling gas shell is physically more re- 
alistic, and it certainly does not offer any advan- 
tages over the infalling gas column that we have 
chosen for most of our infall experiments. 

The formation of rings in the collapse of the 
spherical gas shells is not a big surprise. Ring 
formation in the collapse of rotating spherical 
systems has been documented by several work- 
ers, with transient rings forming in dissipation- 
less collapses ( Miller &: Smith 1979| ), and more 
stable rings forming in the collapse of rotat- 
ing isothermal gas clouds ( p?ohline 198C| ; Boss & 
Haber 1982| ) . The equilibrium of self-gravitating 



isothermal gas rings has been investigated from 
a theoretical perspective as well ( Pstriker 1964 ). 
Unlike the rings seen in dissipationless collapses 
([Miller fc Smith 1979| ), the rings we see are not 
transient structures that are formed en route to a 
more permanent morphology. Furthermore, they 
are true rings in that they are not due to or- 
bit crowding, and the particles in them are not 
constantly moving in and out of the ring con- 
figuration. Our rings are consistent with those 
seen and studied in isothermal gas cloud collapse 
scenarios, and are self-gravitating: the mass per 
unit length is much greater than the equilibrium 
mass per unit length for isothermal rings (Ps-| 
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triker 1964, eq. 156). For G7, this quantity is 
~ 4 times, and for G8, it is ~ 7 times the equi- 
Hbrium value at the last time-step, which cor- 
responds to ~ lltff for G7 and ~ 8ts for 08 
(where is the free-fall time of the initial spher- 
ical shell). However, collapse is most likely inhib- 
ited because our numerical resolution (as set by 
the softening length, which is 0.5 kpc) is of the 
order of the ring width. Our values of a and /?, 
the ratios of the thermal and rotational energies 



to the gravitational energy respectively ( poss fc 
Ha her 1982| ), are 0.02 and 0.12 for 07, and 0.02 
and 0.07 for 08. The value of a is well below 
the maximum value that will produce a collapse. 
For 07, the /? value falls within, whereas for 08, 
the P falls below the 0.1 — 0.5 range seen by Boss 
& Haber (1982) for isothermal rings. However, 
a lower /? limit for isothermal rings has not been 
established as yet, and it should be remembered 
that our initial configuration is a spherical shell 
rather than a cloud. 

4. Discussion 

4.1. SPH vs. Sticky Particles 

Our primary aim with sticky particle simu- 
lations (TR and TRJB) was to determine the 
range of input parameters that worked, and to 
investigate the effect of the infall/merger on the 
primary disk. We succeeded in achieving these 
objectives, and we did not put much faith in the 
sticky particle results where the structure of the 
counterrotating disk was concerned. The SPH 
simulations presented here confirm that we were 
justified in our mistrust of the sticky particle re- 
sults. Shocks and viscosity paint quite a different 
picture of counterrotating disk structure. 

The most visible differences in the disks formed 
with SPH are their thinness and radial structure. 
Sizes are comparable to the sticky particle disks, 
although the SPH disks do on average tend to 
be a little smaller. However, they are all signif- 
icantly thinner than their sticky particle coun- 
terparts. We found the sticky particle disks to 
be on average thicker than the primary disks. 
With SPH the reverse is true. The velocity field 
diagrams also show that all the counterrotating 



disks are considerably colder than the primary 
disks. This was true even with the sticky parti- 
cle disks. 

There is more structure visible in the SPH 
disks. A comparison of the final panels in Figs. 
|, and I with Figs. 1, 4, 6 (except 7) of TR 
and Figs. 2-4 and 8 of TRJB illustrates the dif- 
ferences in the nature of the SPH disks. None 
of the disks produced with sticky particles had 
radial profiles approaching an exponential mass 
distribution. We were able to come close with 
SPH (runs 02 and 04). The counterrotating disk 
produced in 04 shows clear evidence of lopsided- 
ness (Fig. ^), and warps are seen in nearly all of 
the thin disks formed with SPH. 

The time-scale of the counterrotating disk 
formation is significantly smaller (by ^ 30%) 
with SPH. A comparison of simulations 01-04 
with their sticky particle counterparts proves 
this. The sticky particle run corresponding to 
01 produced a counterrotating disk in ^ 6 Gyr 
(TR), whereas 01 and 02 take ;^ 4 Oyr. For 
NOC 4138, sticky particle simulations were able 
to produce a counterrotating disk in ;^ 3 Oyr 
(TRJB), whereas 03 and 04 achieve the same 
result in ^ 2 Oyr. Most of the decrease in the 
time required is due to the shorter columns of 
gas used in the SPH simulations, but a compar- 
ison of the evolution of the sticky particle (we 
ran a sticky particle simulation with the same 
inputs for comparison) and SPH gas also shows 
that the latter loses its kinetic energy and angu- 
lar momentum more rapidly. We estimate that 
the difference in the time-scales due to the gas 
dynamics alone is of the order of 5-10%. The in- 
falling gas is captured and assimilated faster into 
the primary disk because of the more efficient 
viscous dissipation in SPH. 

Strangely enough, the infalling gas does not 
collapse and fragment like it did with sticky par- 
ticles. The basic character of dissipation ap- 
pears to be different: in sticky particle dissipa- 
tion, head-on collisions were not quite so dissipa- 
tive, and there was a comparable amount of dis- 
sipation regardless of how much collisional sup- 
port the gas had; in SPH, whenever collisions be- 
tween gas particles cause shocks, the dissipation 
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is drastic and the gas quickly loses most of its 
collisional kinetic energy. As a result, as long as 
the SPH gas motion does not involve any head- 
on collisions between opposing gas streams, the 
gas dissipates very little energy. The tendency 
of SPH to cause rapid loss of kinetic energy and 
angular momentum in counterstreaming gas has 



been noted by other users of SPH as well (Mihos 
19^0), and it is possible that SPH is too dissipa- 
tive in nature (the artificial viscosity too high) 
under these circumstances. 

4.2. Source of Infalling Material 

While it is fine to theorize about the forma- 
tion of counterrotating disks by gas infall or even 
gas-rich dwarf infall, the evidence for such mate- 
rial in the intergalactic medium is still quite thin, 
which is the main reason why there are no real- 
istic models of infall available. The possibility 
of undetected molecular gas in the IGM around 
spirals provides one glimmer of hope. There 
has been more evidence recently for a significant 
fraction of the neutral ISM in our galaxy being 
i?2-bearing molecular clouds ( pixon, Hurwitz~fc 
Bowyer 1998| ). The general case for ongoing gas 
infall in spiral galaxies is also getting stronger, 
with recent evidence from IRAS data (albeit still 
somewhat controversial), that massive star for- 
mation rates are independent of spiral type and 
that the median gas recycling time is ~ 1 x 10^ 
yr ( pevereux &: Hameed 1997| ) , pointing to a sus- 
tainable outside source of gas. 

"Galaxy harassment" has been recently sug- 
gested as a means of spiral evolution in clusters 
( pVIoore et al. 1996 ). Galaxy harassment refers 
to multiple, high-speed encounters that trans- 
form small spirals into dE/dSph, but more im- 
portantly, leave giant tidal debris arcs and tails 
that could provide fuel for quasars as well as the 
raw material for producing counterrotating disks. 
The recent discovery of a giant debris arc in the 
Coma cluster ( Trentham &: Mobasher 1997] ) sup- 



ports this theory. Harassment of the debris tails 
themselves would create tidal shocks leading to 
formation of dwarf galaxies in the tidal tails, as 



dence of counterrotating disks in clusters should 
be higher than in the field, something that can 
be tested for in the future as more systematic 
surveys of counterrotation in clusters (especially 
the more distant ones) are undertaken. 

The main difficulty posed by dwarf mergers is 
the heating experienced by the primary disk due 
to the large amounts of dark matter currently be- 
lieved to be associated with most dwarf galaxies. 
There is still some doubt whether the evidence 
for the dark matter in dwarf galaxies is trustwor- 
thy or not, and if vindicated, dwarf spheroidal 
satellites without dark matter ( |Kroupa 1997 ) 
may provide a better chance of obtaining coun- 
terrotating disks from dwarf mergers. The in- 
stances of counterrotating bulges in spiral galax- 
ies (e.g. NGC 7331, prada et al. 1996|) and 



counterrotating cores in ellipticals certainly point 
to merger events (iBalcells fc Quinn 199q) . We 
are able to create a counterrotating bulge with 
our gas-rich dwarf model when the dwarf is suffi- 
ciently dense and massive (simulations SI and S2 
in TR), but the dissipationless (dark and stellar) 
matter in the dwarf is a severe liability to the 
survival of the primary disk. 

4.3. Gravitational Influence of Various Com- 
ponents 

To compare the gravitational influence of the 
three components (halo, disk and gas) in the disk 
plane, we have plotted the gravitational force ex- 
perienced by a test particle of unit mass in the 
mid-plane of the disk (z = 0) as a function of 
radius in Fig. 21. This is shown at the begin- 



shown also in numerical simulations ( Barnes fc 
Hernquist 1992| ). This would mean that the inci- 



ning and end of the simulation, for simulations 
G2 (cold primary) and G4 (hot primary). The 
gravitational forces experienced by the test par- 
ticle due to the halo particles, the primary disk 
particles, and the gas particles, are shown sepa- 
rately. The forces are further separated into their 
azimuthally averaged radial (in plane of disk) and 
vertical components. Values for the outer halves 
of the disks are less reliable due to low particle 
numbers, and for G2 the calculation of the force 
components is further complicated by the sub- 
stantial inclination of the primary and gas disks 
at the end of the simulation. We have attempted 
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to correct for this inclination as best we can by 
first applying a coordinate transformation to re- 
set the disk plane to the xy-plane before doing 
the force calculations. 

For the cold primary (G2), initially {t = 1.5 
Gyr) the primary disk exerts the dominant ra- 
dial force in the inner half of the disk system. 
The disk and halo radial forces are almost equal 
in the outer half. The vertical forces due to all 
components are initially an order of magnitude 
lower than the radial forces. As expected, the in- 
fluence of the infalling gas is negligible initially. 
This is because the gas is several disk radii away 
at the beginning of each simulation. At the end 
of the simulation (t = 6.0 Gyr), the halo-disk 
force distribution with radius does not change 
significantly, but the magnitude of the forces is 
lower on average by ^ 7%. The radial force due 
to the gas is now comparable to the halo and disk 
radial forces. In the first couple of kpc, the gas 
radial force dominates. The vertical influence of 
the halo has increased and that of the primary 
disk has diminished by a considerable amount 
(~ 50 — 60%), especially in the outer half of the 
disk plane. The gas now exerts a vertical force 
that is dominant in the inner half of the disk, and 
is intermediate between the halo and disk verti- 
cal forces in the outer half. On average, though, 
the vertical forces are still an order of magnitude 
lower than the radial forces. 

For the hot primary (G4), the initial {t = 0.5 
Gyr) breakdowns are similar although the halo is 
even less dominant in the inner half of the disk 
plane due to its high core radius (5 kpc). The 
vertical forces are as before an order of magni- 
tude lower than the radial forces, and the gas 
influence is negligible initially. At the end of the 
simulation (t = 2.5 Gyr), the influence of the gas 
is close to that of the halo in the inner regions, 
but is signiflcantly less than the disk and halo 
forces in the outer disk. The vertical influence 
of the halo remains much the same, whereas the 
vertical force due to the primary disk increases 
in the inner half. The gas vertical force rivals 
that due to the halo throughout, and is not the 
dominant vertical influence as was the case for 
the cold primary. The gas mass fraction (com- 



pared to the primary disk mass) is much less in 
the hot primary, and so it is not surprising that 
the gas has a smaller influence both radially and 
vertically. 

4.4. The Nature of the Gas Disks 

The final rotation curves and mid-plane vol- 
ume densities of the counterrotating gas disks at 
the end of simulations G2 and G4 are shown in 
Fig. 22. The rotation velocity and velocity dis- 



persion values (azimuthally averaged) are shown 
in the plots at the top, with the SPH-computed 
total mid-plane gas density (azimuthally aver- 
aged) at each radial location is shown with filled 
circles in the bottom plots. This mid-plane den- 
sity is compared to /jqo, the unperturbed central 
density for the plane-parallel obtained by adding 
up the contributions due to the gas pressure and 
the (extern al) gravitational pressure as pe r equa- 
tion (9) in [Elmegreen fc Elmegreen 197£ : 



Poo = (fcxt + ^TTGa^)/c^. 

This "expected" mid-plane density is plotted as 
the open circles, and the second term on the right 
of the above equation, due to the gas surface den- 
sity, is plotted as the open squares. The gravita- 
tional pressure in the first term is obtained from 
the combined halo, primary disk and gas vertical 
forces shown in Fig. ^ 

For the cold primary (G2) gas disk, the ro- 
tation velocity is steep initially and levels off at 
i? ~ 3 kpc, whereas the hot primary (G4) gas 
disk's rotation curve has a more gentle slope to- 
wards its flat portion. This is due to the higher 
halo core radius for the NGC 4138 model, which 
reduces the halo's influence in the inner few kpc. 
The heights of the flat portions of the rotation 
curves are comparable for the gas disks in G2 and 
G4, but the G2 gas disk clearly shows a higher 
average velocity dispersion. The counterrotat- 
ing disk formed in the cold primary is evidently 
somewhat hotter than the one formed in the hot 
primary. 

The value of poo i the expected mid-plane den- 
sity obtained by adding up the two contributions 
from the gravitational and gas pressures, com- 
pares well with the actual computed mid-plane 
gas density, although there are signiflcant differ- 
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ences. For the cold primary, the expected value 
is consistently less than the actual density except 
for the outer disk. In the case of the hot primary, 
the reverse is true. The differences seen in the 
computed and expected densities are probably 
due to errors in the surface density and gravita- 
tional pressure values obtained from the particle 
data. Sources of uncertainties include the cutoff 
height assumed for the gas disk, the correction 
applied for the disk center-of-mass and inclina- 
tion (especially in G2), and lower particle num- 
bers at higher radii. The curves for the gas pres- 
sure (surface density) term indicate that the gas 
pressure is the dominant contributor for the gas 
density for the hot primary, whereas for the cold 
primary the gravitational pressure is the domi- 
nant influence, especially in the outer disk where 
the gas surface density falls off considerably. 

4.5. Heating of Primary Disk 

The initial and final rotation curves for the 
primary disks after gas infall are shown in Fig. 
for both the cold (G1-G2) and hot (G3-G4) 
primary simulations. The final rotation veloci- 
ties and velocity dispersions are consistent with 



those observed for SO/Sa galaxies (e.g., Seifert fc 
Bcorza 1996; Rix et al. 1992| ) since frot/c ^ 2 



rather than ^ 1 as is true for Sb and later types. 
Furthermore, it is clear from Fig. ^ that both the 
cold and hot primaries end up looking like SOs 
after the formation of the counterrotating disk. 
The velocity dispersion ranges from ~ 100 to 
~ 150 km/s for the final disks in either case, with 
the rotation velocities dropping by ~ 50 — 100 
km/s for the cold primary, and ~ 25 — 30 km/s 
for the hot primary. In other words, the cold pri- 
mary fares much worse than the primary that was 
hot (SO/Sa) to begin with. In the same vein, the 
outer portions of the primaries, which are colder 
initially (since the velocity dispersion is initially 
proportional to the surface density), are heated 
up more than the inner regions. 

As mentioned above, the heating seen in simu- 
lations G1-G4 is more than was seen for the cor- 
responding sticky particle simulations, because 
the infalling gas columns are shorter in the SPH 
simulations. However, the fact that colder disks 



are more susceptible to heating is independent of 
the actual amount of heating experienced. Al- 
though this vulnerability of cold disks is hardly 
a surprise, it does strengthen the argument for 
predicting that most instances of massive coun- 
terrotating disks are likely to be found in SO/Sa 
galaxies. This is consistent with the observations 
to date. 

4.6. Tilting of Disks 

It is apparent from most of the simulations 
presented above that the resulting counterrotat- 
ing gas disk is significantly tilted with respect to 
the initial symmetry plane of the primary disk 
(the z = plane). This is not a phenomenon 
new to our SPH simulations. It was observed 
in the sticky particle simulations as well and we 
have commented on it previously (TRJB). It is 
worth reiterating here that the tilting experi- 
enced by the gas disk is in phase with the tilting 
of the primary stellar disk, occurs over relatively 
short time-scales (order of dynamical time of pri- 
mary disk) compared with the overall time-scale 
of each simulation, and is in response to the com- 
bined torques exerted by the halo and the gas. 

Even though the net gravitational torques and 
the angular momentum vectors are initially all 
aligned with one another, the forces experienced 
by individual particles subsequently are by no 
means all in the initial plane of the primary disk. 
This is especially true for the gas particles, which 
undergo gas dynamical interactions. For instance 
in G2, which shows the maximum final disk in- 
clination, the side views for t = 3.0 and t = 4.0 
Gyr (Fig. ^) show an excess of gas particles below 
the disk plane (i.e., in the lower half of the pan- 
els, which are centered on z = 0) on the left side 
of each panel. This excess mass concentration 
below the plane on one side is consistent with 
the subsequent tilt of the gas disk (which tilts in 
phase with the primary disk) seen in later panels. 

Comparisons of the total torques exerted by 
the halo and the infalling gas on the primary disk 
with the rate of change of the angular momen- 
tum vector (TRJB) of the disk show that the two 
quantities are consistent with each other. The 
change in the value of the z-component of the 
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total angular momentum vector, which is by far 
the dominating component, is less than 1% over 
the course of the simulation in most cases, and 
up to ~ 3% in the extreme cases. 

4.7. Star Formation 

Our SPH code does not incorporate star for- 
mation as yet. Previous studies suggest that in 
interactions of galaxies where substantial amounts 
of gas are involved, the star formation rates are 
very modest for more than 90% of the interac- 
tion history, with almost all of the star formation 
occurring in a rapid starburst after the gas has 
dissipated most of its kinetic energy and formed 
dense regions in the centers of the galaxies (Mi 



tios, Richstone & Bothun 1992|; 


Mihos, Bothun 


fe Richstone 1993; 


Mihos k Hernquist 1994a, 



1994b). The pre-starburst phase of the star for- 
mation in these simulations had no discernible 
effect on the final structure of the gas, since the 
depletion rates were very low. In our case, the 
longer time scales involved in our retrograde sim- 
ulations may deplete significantly more gas over 
the course of the simulation, but other than re- 
ducing the intensity of the central starburst, this 
is unlikely to have a profound effect on the struc- 
ture of the gas disks formed as long as there is 
little or no prograde gas in the primary disk. If 
there is a considerable amount of preexisting pro- 
grade gas, then it is possible that star formation 
resulting from the head-on collision of the coun- 
terstreaming gas particles will prevent a signifi- 
cant amount of gas from falling to the center due 
to rapid conversion to stars. This may alleviate 
the problem of excessive dissipation that causes 
the gas to collapse to the center of the primary 
in our simulations that include prograde gas (G5 
and G6). The presence of star formation may 
also be instrumental in producing counterrotat- 
ing radial profiles that are closer to the primary 
radial profiles (i.e., closer to exponential profiles). 
We plan to include star formation in our code in 
the future, but for now we have excluded it to 
save precious CPU time. 



5. Conclusions 

There are some notable differences between 
the characteristics of the counterrotating disks 
resulting from SPH and those obtained with 
sticky particle simulations by TR and TRJB. The 
SPH disks are very thin compared to their sticky 
particle counterparts and compared to the pri- 
mary disks. They also show evidence of spiral 
structure, and their size and radial mass distribu- 
tion are quite sensitive to the input parameters, 
particularly those that affect the initial angular 
momentum of the gas. 

Other differences in the SPH results include 
the lack of clumping of the infalling g cLS, Si prob- 
lem that was quite severe with our sticky particle 
simulations, and the shorter time-scales for disk 
formation. 

Although it is easy to produce thin counter- 
rotating disks with gas infall, it is not so easy 
to obtain exponential radial profiles. The ini- 
tial angular momentum of the gas has to be low 
enough, and some combination of other processes 
such as prograde gas, star formation and energy 
feedback from massive stars may be necessary to 
produce counterrotating disks that are very sim- 
ilar to the primary disks. Currently there is no 
evidence to indicate that counterrotating disks 
have predominantly exponential profiles, so this 
is not necessarily a problem. 

In general, the process that dumps a massive 
counterrotating disk in a cold primary spiral, es- 
pecially if it is a minor gas-rich merger but even 
if it is gas infall that occurs over a few dynamical 
times (^ lOtdyn)) is likely to heat up the primary 
substantially and change its type to an SO/Sa 
galaxy. On the other hand, if the primary is al- 
ready an SO/Sa galaxy to begin with, then it can 
acquire a counterrotating disk without changing 
its type significantly. The fact that most of the 
currently known instances of massive spiral coun- 
terrotating disks are in SO/Sa galaxies is there- 
fore a selection effect rather than an accident. 

The presence of primordial prograde gas in the 
primary has a drastic effect on the retrograde gas 
that comes in contact with it. Neutralization of 
the angular momenta is rapid, with both the pro- 
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grade and retrograde gas particles ending up in 
the center of the primary within a few dynamical 
times. This may be an indication of a problem 
with SPH that causes over-dissipation in coun- 
terstreaming gas flows, at least in the absence of 
star formation. The inclusion of star formation 
and energy feedback from super novae will most 
likely yield significantly different results in such 
situations. 

A retrograde-rotating, infalling gas shell pro- 
duces a counterrotating bulge and flat outer ring, 
but is unable to produce a counterrotating disk in 
the proper sense. The size of the ring is well cor- 
related with the angular momentum of the shell. 
The formation of the ring is consistent with pre- 
vious studies of collapsing isothermal gas clouds. 
These studies also suggest that significantly hot- 
ter gas with lower angular momentum is neces- 
sary to produce a counterrotating disk with this 
model. 

We hope to test our results further in the near 
future with the addition of thermal effects and 
star formation to our SPH code. 
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Fig. 1. — Run Gl - counterrotating disk formation by retrograde infall of gas onto a cold primary disk. 
The top (upper 8 panels) and side (lower 8 panels) views are shown. The primary disk is shown only 
in the first panel of each view for clarity. The first two panels for each view are 160 kpc wide, with 
subsequent panels having the zoom factors indicated. Time is in Gyr. 

Fig. 2. — Run G2 - infall of gas with lower angular momentum than in Gl. Panel widths and time units 
are the same as in Fig. |l|. The radial distribution of the gas is much closer to an exponential disk rather 
than a ring as in Gl. 

Fig. 3. — A comparison of the primary disk for infall simulation G2 at t = 1.5 Gyr (top row) and t = 6.0 
Gyr (middle row) and the counterrotating gas disk at t = 6.0 Gyr (bottom row), with the top views on 
the left and the side views on the right. The width of each panel is 40 kpc. 

Fig. 4. — Run G2 - velocity fields of the primary disk (left) and counterrotating gas (right). Panel width 
is 40 kpc, and a fraction of the particles is sampled for clarity, with sampling proportional to radius. 

Fig. 5. — A comparison of the radial density profiles of the counterrotating disks for runs Gl (dotted 
line) and G2 (dashed line), measured at the final timesteps. The surface particle number-density is 
calculated by counting and azimuthally averaging the total number of particles in columns of unit area 
and half-height (1 kpc^ and 1 kpc respectively), and plotted as a function of radius R. A correction is 
applied for the disk inclination if necessary. The radial profile of the primary disk at the end of G2 (the 
Gl primary profile is very similar) is also plotted for reference (solid line). 

Fig. 6. — A comparison of the primary disk thickness for runs Gl and G2. The mean half-thickness z is 
plotted as a function of radius R. The solid curve shows the initial primary disk at t = 1.5 Gyr, and the 
isolated disk curve shows the thickness after evolving the primary disk without any gas accretion for the 
same time period. 

Fig. 7. — Run G3 - counterrotating disk formation by retrograde infall of gas for NGC 4138 primary 
disk. The top (upper 8 panels) and side (lower 8 panels) views are shown. The primary disk is shown 
only in the first panel of each view for clarity. The first two panels for each view are 80 kpc wide, with 
subsequent panels having the zoom factors indicated. Time is in Gyr. 

Fig. 8. — Run G3 - velocity fields of the primary disk (left) and counterrotating gas (right). Panel width 
is 17 kpc, and a fraction of the particles is sampled for clarity, with sampling proportional to radius. 

Fig. 9. — Run G4 - infall of gas with lower angular momentum and a shorter column of gas than in G3. 
Panel widths and time units are the same as in Fig. |^ The radial distribution of the gas is much closer 
to an exponential disk rather than a ring as in G3. 

Fig. 10. — A comparison of the radial density profiles of the counterrotating disks for runs G3 (dotted 
line) and G4 (dashed line), measured at the final timesteps. The surface particle number-density A'^ is 
calculated as in Fig. |5| and plotted as a function of radius R. The radial profile of the primary disk at 
the end of G3 (the G4 primary profile is very similar) is also plotted for reference (solid line). 

Fig. 11. — Run Dl - gas-rich dwarf merger M2 from TRJB repeated with SPH. Panel widths and time 
units are the same as in Fig. 0. 
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Fig. 12. — Run Dl - velocity fields of the primary disk (left) and counterrotating gas (right). Panel width 
is 20 kpc, and a fraction of the particles is sampled for clarity, with sampling proportional to radius. 



Fig. 13. — Run D2 - partial SPH run of gas-rich dwarf merger M2 from TRJB. Panels are 10 kpc wide. 



Fig. 14. — Run G5 - A prograde ring of gas is included in the primary disk. Both the prograde and 
retrograde gas particles are shown. The primary stellar disk is not shown. Panel widths and time units 
are the same as in Fig. |9|. 



Fig. 15. — Run G5 - Velocity fields of the primary prograde gas (left) and incoming, initially retrograde 
gas (right) shown for the first four timesteps of Fig. 14. Each panel is 4 kpc wide. Only a fraction of the 
particles is shown for clarity. 



Fig. 16. — Run G6 - A prograde disk of gas is included in the primary disk with all other input parameters 
the same as in G4. Both the prograde and retrograde gas particles are shown. The primary stellar disk 
is not shown. Panel widths and time units are the same as in Fig. 0. 



Fig. 17. — Run G6 - velocity fields of the corotating (left) and counterrotating (right) gas at t = 1.1 
(top), t = 1.5 (middle) and t = 1.7 (bottom) Gyr. The panels on the left are 1.5 kpc wide, and the panels 
on the right are 4 kpc wide. Only a fraction of the gas particles are shown for clarity. 



Fig. 18. — Run G7 - Retrograde infalling (initially spherical) gas shell. The top views (looking along the 
z-axis) are shown on the left and the side views (looking along the y-axis) are shown on the right. The 
shell is initially 10 kpc thick and its initial radius is 30 kpc. The primary stellar disk is not shown. Panel 
widths for the first two panels in each view are 80 kpc, with subsequent panels having the zoom factors 
indicated. 
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Fig. 19. — Run G8 - Retrograde infalling shell with half the initial angular momentum compared to the 
shell in G7. The top views are shown on the left and side views on the right. The primary stellar disk is 
not shown. Panel widths for the first two panels in each view are 80 kpc, with subsequent panels having 
the zoom factors indicated. 



Fig. 20. — Top and side views of the counterrotating disk formed in G4 showing evidence of lopsidedness 
(top) and warps (bottom). 



Fig. 21. — The azimuthally averaged gravitational force due to the halo (circles), disk (squares) and gas 
(triangles) particles, at points located in the plane of the disk at the plotted radii, with the component 
in the plane shown as the solid line, and the component perpendicular to the plane shown as the dotted 
line. The plots are shown for the initial (top) and final (bottom) configurations, with results for the cold 
primary simulation G2 being shown on the left and the hot primary simulation G4 shown on the right. 
The force values are in simulation units, time in Gyr. 



Fig. 22. — The final rotation velocity (open circles) and velocity dispersion (open squares) plots (top), 
and the mid-plane gas density (bottom) for the gas disks formed at the end of simulations G2 (left) and 
G4 (right). In the gas density plots, the actual average mid-plane gas density computed using SPH is 
plotted as the filled circles, whereas poo, the expected mid-plane density (see text), is shown as the open 
circles. The contribution to the gas density due to the gas pressure (see text) is plotted as the open 
squares. The density values are in simulation units. 



Fig. 23. — A comparison of the rotation velocity v^ot and velocity dispersion cr as a function of cylindrical 
radius R of the initial primary disks (top) with the final primary disks (middle and bottom) for gas infall 
simulations with a cold primary (left) and hot primary (right). The open circles show the mean rotation 
velocity and the open squares the velocity dispersion. The velocity unit is ~ 1000 km/s. Both velocities 
for the final disks have been corrected for disk inclination, but values for R ^ OAR^, where Rd is the 
maximum value of R plotted, are not reliable. The filled triangles show the ratio Vrot/c (limits shown on 
the right axis) as a function of radius. 
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